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ABSTRACT 

We report results of iV-body simulations of isolated star clusters, performed up to the 
point where the clusters are nearly completely dissolved. Our main focus is on the post- 
collapse evolution of these clusters. We find that after core collapse, isolated clusters 
evolve along nearly a single sequence of models whose properties are independent of 
the initial density profile and particle number. Due to the slower expansion of high-iV 
clusters, relaxation times become almost independent of the particle number after 
several core collapse times, at least for the particle range of our study. As a result, the 
dissolution times of isolated clusters exhibit a surprisingly weak dependence on N. 

We find that most stars escape due to encounters between single stars inside the 
half-mass radius of the cluster. Encounters with binaries take place mostly in the 
cluster core and account for roughly 15% of all escapers. Encounters between single 
stars at intermediate radii are also responsible for the build up of a radial anisotropic 
velocity distribution in the halo. For clusters undergoing core oscillations, escape due 
to binary stars is efficient only when the cluster center is in a contracted phase. Our 
simulations show that it takes about 10 5 iV-body time units until the global anisotropy 
reaches its maximum value. The anisotropy increases with particle number and it seems 
conceivable that isolated star clusters become vulnerable to radial orbit instabilities 
for large enough N. However, no indication for the onset of such instabilities was seen 
in our runs. 

Key words: methods: numerical - celestial mechanics, stellar dynamics - globular 
clusters: general. 
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1 INTRODUCTION 

It is one of the great challenges of numerical astrophysics 
to understand the dynamical evolution of globular clusters 
and other collisional systems (as e.g. galactic nuclei; galaxy 
clusters). So far however, it is not possible to perform simu- 
lations of such rich systems by direct summation techniques 
and one has to rely on a statistical modeling. This is still 
the case, despite recent progress in the development of faster 
computers, like e.g. the GRAPE computer series (Makino 
2001, 2002). 

Statistical modeling of globular clusters includes differ- 
ent techniques like Fokker-Planck, Monte Carlo and Gaseous 
model simulations. There has been substantial progress in 
recent years in improving the reality of such models and 
comparing them with iV-body models or comparisons be- 
tween different statistical methods show good agreement in 
many cases (Takahashi & Portegies-Zwart 2000, Giersz & 
Spurzem 2000). Nevertheless the validity of such codes for 
real globular clusters should be checked further, due to the 
complexity of the physical processes involved and the fact 



that the gap between the largest published A^-body models 
(which are of order 3 • 10 4 ), and the number of stars in real 
globular clusters (roughly N — 10 6 ) is still large. 

Isolated star clusters are an ideal environment to test 
star cluster evolution. These systems are sufficiently simple, 
but nevertheless show many phenomena that are also seen 
in real globular clusters. Although simulations of isolated 
clusters are not directly applicable to real stellar systems, 
they are nevertheless highly relevant. Borrowing an exam- 
ple from stellar evolution, such calculations play the role 
that polytropic models have played there historically. While 
rather unrealistic physically, they did shed considerable light 
on the basic aspect of stellar structure, and in addition they 
served as useful comparison material when trying to inter- 
pret the results of the far more detailed numerical evolution 
models. Models for isolated clusters have the same dual role 
in stellar dynamics, they show some of the major physical 
effects occurring in long-term dynamical evolution, and pro- 
vide templates against which more detailed models can be 
interpreted. A large amount of research has therefore been 
devoted to the study of isolated clusters. 
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It has been known for a long time that star clusters un- 
dergo a contraction of their center as a result of heat transfer 
from the cluster core to the halo (Antonov 1962, Lynden-Bell 
& Wood 1967). For isolated, single-mass clusters, this core 
collapse takes roughly 17 initial relaxation times until com- 
pletion (Takahashi 1995), and the core collapse evolution 
is thought not to change very much in the presence of an 
external tidal field (Giersz & Heggie 1997). Lynden-Bell & 
Eggleton (1980) found that the core collapse proceeds from 
the outer to the inner parts and leaves behind a power-law 
debris of material with density exponent in the range a — — 2 
to -2.5. Although core collapse strongly influences the struc- 
ture of the inner parts, one usually assumes that the overall 
properties of the system are determined more by the half- 
mass radius, and the core adjusts itself in post-collapse so 
as to balance the tendency of the halo to recollapse. 

Core collapse will be reversed if a central energy source 
is present (Henon 1975), and it is normally assumed that 
binary stars, which are either primordial or form during 
core collapse and harden as a result of encounters with field 
stars, provide this heat source (Aarseth 1971, Heggie 1975) 
in real star clusters. The number of binaries necessary to 
drive the cluster evolution is found to be quite small. Good- 
man (1984) estimated from his gaseous model calculations 
that their number drops with the number of cluster stars 
as JV;, oc iV -0 ' 3 and that in the mean only Nt = 0.5 bina- 
ries are necessary to drive the evolution of a cluster with 
N — 10 6 stars. This has been partly confirmed by Giersz & 
Heggie (1994b) in their iV-body simulations of small- N clus- 
ters. They found that only few binaries are present in the 
core, but their number seemed to increase with the parti- 
cle number. They noted however that some of their binaries 
might be temporary ones which contribute nothing to the 
energy generation of the core. 

Stars scattered out of the cluster center create an radi- 
ally anisotropic velocity profile in the cluster halo (Larson 
1970). Cohn (1984) found that a cluster starts to become 
anisotropic in its outer parts already early on in the pre- 
collapse phase, while Spitzer & Shapiro (1972) have argued 
that the anisotropy extends down into the core as a conse- 
quence of its collapse. This was later confirmed by Takahashi 
(1995) and Drukier et al. (1999) for the late stages of core 
collapse. For post-collapse clusters, Takahashi (1996) found 
that they are isotropic in their center and the anisotropy 
increases monotonically towards the halo. 

Encounters of stars in the cluster center also lead to 
the escape of stars. Henon (1960) has shown that stars in 
isolated systems do not escape by distant encounters, but 
only due to close encounters with other stars. Despite their 
small numbers, it is generally believed that binaries play an 
important role in the production of escapers and contribute 
to the development of an anisotropic velocity profile in the 
cluster halo. 

So far, most simulations of isolated star clusters were 
concerned only with the early stages up to core collapse or 
the immediate re-expansion phase. Nearly nothing is known 
about how clusters evolve when substantial mass-loss sets in. 
One reason for this neglect is that the escape of stars from 
isolated clusters happens very slowly and on a timescale 
which increases rapidly as the clusters expand. In addition, 
large- angle scatterings are neglected in Fokker-Planck meth- 



ods and can only be studied by iV-body or Monte Carlo 
simulations. 

Hence, although highly desirable, it is difficult to ver- 
ify the predictions of Fokker-Planck or Gaseous models by 
direct iV-body simulations. We therefore present in this pa- 
per the first iV-body simulations of medium sized isolated 
star clusters which follow the evolution until near complete 
disintegration. 



2 DESCRIPTION OF THE RUNS 

We simulated the evolution of clusters containing between 
TV = 128 and 8192 equal-mass stars, increasing the parti- 
cle numbers by successive factors of 2. Our clusters followed 
Plummer profiles initially and we did not impose a cut-off 
radius when the Plummer models were created. Small-iV 
runs were made with the collisional Aarseth iV-body code 
NBODY6 (Makino & Aarseth 1992, Aarseth 1999) on a 
Pentium PC workstation. The 8K run was made with the 
NBODY6++ code (Spurzem 1999, Spurzem & Baumgardt 
2002), which is a parallelized version of NBODY6. It was 
made on a Cray T3E parallel computer using 8 processors. 
Though these are not large by the standards of current 
more realistic simulations, the computations are at least as 
time-consuming. 

During the runs, we recorded the positions and veloci- 
ties of all stars at regular intervals. In the first 100 TV-body 
time units, the output time 8T was equal to one initial cross- 
ing time. Later, time-intervals were equally spaced in log T. 
At each data output, we determined the number of bound 
stars, the lagrangian radii of the members (measured rela- 
tive to the density center) and the velocity structure of the 
bound stars. We also calculated a scaled time according to 

?t nbody f rh (0)\ 3/2 ,. m 

where ru is the half-mass radius of the bound cluster. This 
time takes the expansion of the cluster into account (but 
not the mass-loss) , and measures how much time has passed 
in comoving coordinates. Data was also stored every 2.828 
timeunits in scaled time. 

The runs were performed for a maximum time of 
Tmux = 10 16 TV-body units. In order to follow the evolu- 
tion up to the point where nearly all stars are unbound, 
we had to rescale the whole system when its size became 
too large. Otherwise the simulations would have come to a 
halt due to numerical difficulties in NBODY6. Rescaling was 
done by calculating the sum of the potential energies of all 
stars, excluding the internal energy of regularised binaries. 
The positions and velocities of all stars were then rescaled 
such that the total potential energy after rescaling was equal 
to Ep t — ^1/2. The cluster centers were also moved back 
to the origin and the mean cluster velocity was subtracted 
from all stars. This recentering did not influence the clus- 
ter membership. In general, the whole integration could be 
covered by two rescalings. 

Bound stars were determined in the following way: We 
first calculated the potential energy of each star with respect 
to all other stars. Kinetic energies were then calculated rel- 
ative to the motion of the cluster center, which was taken 
from a previous iteration, or, at the start of the run, assumed 
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Table 1. Details of the performed TV-body runs. 



N 


Nsim 


TFinh 


< NTmax > 


Tec 


THalf 


128 


64 


2.8 -10 12 




70.7 


6.75 ■ 10 4 


256 


64 


1.3 10 14 




104.7 


6.70 ■ 10 4 


512 


40 


1.2 10 15 


12.7 


169.7 


7.56 ■ 10 4 


1024 


20 


1.0 -10 16 


28.5 


360.3 


9.05 ■ 10 4 


2048 


10 


1.0 -10 16 


49.4 


587.5 


1.39 ■ 10 5 


4096 


5 


1.0 10 16 


100.5 


1031.2 


2.29 -10 5 


8192 


1 


1.0 -10 16 


182.0 


1839.2 


3.79 -10 5 



to be zero. The mean cluster motion was then determined 
from the bound stars (i.e. those with total energies E < 0) 
and in the following iteration potential energies were cal- 
culated with respect to the bound stars only. This method 
was repeated until a stable solution for the cluster motion 
and the number of bound stars was found. From the cluster 
members, we determined the position of the density cen- 
ter and the lagrangian radii, using the method of Casertano 
& Hut (1985). Unbound stars were not removed from the 
simulations, although for the analysis in this paper we will 
restrict ourselves to the bound stars. We note that through- 
out the paper, the term cluster will be used for the system 
of bound stars only. We also note that using unbound stars 
for the calculation of the potential energy or the position of 
the density center would make little difference to our results, 
since they contribute little to the potential energy and also 
have only a small weight in the calculation of the density 
center. 

Table 1 gives an overview of the simulations. It shows 
the number of cluster stars N, the number Nstm of sim- 
ulations performed, the time TFinh when half the simula- 
tions had stopped (explained below), the number of stars 
still bound at 10 16 iV-body time units and the core collapse 
and half-mass times of the models. 

Runs were terminating before Tmox due to either diffi- 
culties in the numerical integration or because our routine 
for the membership determination couldn't find bound stars. 
This happened for example when a cluster consisted of only 
two or three stars and these were merged by NBODY6 to 
form an hierarchical system. Runs which finish before Tuax 
might bias our results, so results obtained for T ^ Tpi n h 
should be treated with caution. However, Table 1 shows that 
only results for small- TV clusters might be biased and these 
also only for the very late stages of the evolution. 

Throughout the paper we will use iV-body units, so the 
constant of gravitation, initial cluster mass and energy are 
given by G = 1, M = 1 and Ec = —0.25 respectively. 



3 RESULTS 

3.1 Pre-collapse evolution 

We will start our analysis by discussing the pre-collapse evo- 
lution of the clusters. Fig. 1 compares the evolution of dif- 
ferent models after the iV-body time units were divided by 
the half-mass relaxation time (Spitzer 1987) 



T RH = 0.138 



3/2 



(2) 



l 

0.98 
0.96 
0.94 




fm VG InjN 

where N is the initial number of cluster stars, m their 
mean mass, r;, the initial half-mass radius of the cluster, 7 a 
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Figure 1. Pre-collapse evolution of clusters containing TV = 1024 
(dashed) and N = 4096 (solid line) stars. The bottom panel shows 
the evolution of the lagrangian radii, the top panel shows the evo- 
lution of the bound mass. Time is scaled by the initial relaxation 
time. The lagrangian radii evolve in the same way in both models 
until shortly before core collapse. 



constant in the Coulomb logarithm (discussed below), and 
G the constant of gravitation. 

The bottom panel shows lagrangian radii containing be- 
tween 3% to 70% of the bound stars, calculated by aver- 
aging the radii from individual runs. The collapse of the 
clusters and the following re-expansion can be clearly seen. 
Core collapse is achieved in approximately 17 initial relax- 
ation times, which agrees well with published results from 
anisotropic Fokker-Planck runs (Takahashi 1995, Drukier et 
al. 1999) . For most of the pre-collapse phase good agreement 
in the evolution of the lagrangian radii for different N can 
be seen. Significant deviations occur only shortly before core 
collapse, when heat input by binaries starts to influence the 
core evolution. Giersz & Spurzem (1994) found that high- TV 
clusters undergo a core collapse which is deeper and hap- 
pens at later times since binaries in low- TV clusters have to 
create less energy to stabilise the cluster core from further 
contraction. Goodman (1987) found in his gaseous model 
calculations that high- TV clusters have a smaller core radius 
also in the post-collapse phase. As can be seen, these results 
agree very well with the results of our TV-body calculations. 
Fig. 1 shows that the core radii of high- TV models remain 
smaller during the post-collapse phase, while their interme- 
diate radii, probably as a result of the stronger energy gen- 
eration in the core, are larger. 

The outer radii are expanding right from the start of 
the simulations. In agreement with Cohn (1984), we find 
that the halos become anisotropic already during the pre- 
collapse phase (Fig. ^), so the pre-collapse expansion is at 
least partially driven by stars ejected from the cluster cores. 
These stars must be ejected by two-body encounters between 
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single stars, as no binaries are present in the clusters before 
core collapse. 

Comparing the evolution of models with different N al- 
lows us to determine the value of the Coulomb logarithm 
used in eq. 2. From a fit of the lagrangian radii, we typi- 
cally obtained values between 7 = 0.05 and 0.2, depending 
on which radii and particle numbers are compared. These 
values are similar to the value of 7 = 0.11 obtained by 
Giersz & Heggie (1994a) from a comparison of N = 500 
and N = 2000 models. Since the number of our simulations 
is relatively small, these differences are probably within the 
statistical error. 

The upper panel of Fig. [l] depicts the mass-loss. The 
mass-loss rate increases monotonically towards core collapse 
and about 3% of the mass is lost in pre-collapse. There is 
a trend that the 4K-model loses mass at a slightly smaller 
rate in the beginning, and this effect can also be seen if we 
compare other models. The reason could be that stars in 
isolated models escape only by very close encounters, for 
which Henon (1960) found that they lead to a scaling of the 
lifetime proportional to N rather than iV/log(jV). 




T [NBODY] 



3.2 Post collapse expansion 

The post-collapse expansion of an isolated cluster is usu- 
ally assumed to be self-similar. In this case, its evolution is 
characterized by the following two equations: 



dM(t) 

dt 
dr h {t) 

dt 



= -ki 



M(t) 



trh 

r h (t) 

trh 



(3) 



If the variation in the Coulomb logarithm of the relax- 
ation time can be neglected, the solution for the mass-loss 
and the half-mass radius is given by: 



M(t) = Mo it/to)-" 
r h (t) = r hQ {t/t Q ) {2+v)/3 



(4) 



where to is the time of core collapse and v a constant. Good- 
man (1984) obtained values of v ranging from 0.0258 to 
0.100 for the evolution of his homologous sphere. Similar- 
values were also found by Giersz & Heggie (1994b) in their 
JV-body simulations of small- N clusters. On the other hand, 
Drukier et al. (1999) found strong deviations from a r oc i 2//3 
scaling for clusters with N ^ 8000 stars in their anisotropic 
Fokker-Planck simulations, especially for the outer radii. 
They attributed these changes to the buildup of a radially 
anisotropic velocity distribution in the cluster halos. 

Fig. |3] shows the evolution of the lagrangian radii as a 
function of iV-body time in our simulations. We have shifted 
the iV-body time units such that all clusters go into core col- 
lapse at the same time as the 4K model. Due to their smaller 
relaxation times, low-JV clusters need less time to reach core 
collapse and also expand quicker after core collapse, which 
can be seen most clearly in the evolution of the outer la- 
grangian radii. 

A strong increase of the outer radii might be expected 
since the halo density in a Plummer model drops off like 
p oc r -5 , while for an isolated cluster with mass-loss one 
expects to find an equilibrium relation p oc r~ 3,5 (Spitzer & 
Shapiro 1972). The build up of this halo will cause the outer 



Figure 2. Ratio of the relaxation times of the IK and 8K models 
at different radii. The 8K model has a smaller relaxation time in 
the core but a larger relaxation time in the halo. Relaxation times 
are approximately the same at the 10% radius, showing that the 
cluster expansion is driven from around this radius. 



radii to expand and might also explain the increase seen by 
Drukier et al. (1999, their Fig. 6). 

The quick expansion of \ow-N clusters cannot be sus- 
tained indefinitely however, since, as a cluster expands, its 
relaxation time increases with it. Hence, the expansion of 
a low-iV cluster cannot go beyond a point where its relax- 
ation time reaches the relaxation time of an high-iV cluster. 
After some time has passed in post-collapse expansion, one 
would therefore expect that high-iV clusters have smaller 
radii which compensate their large particle numbers, and 
that all clusters expand with the same rate. It can be seen 
in Fig. ^ that this is nearly the case. 

Interestingly, from a comparison of different models, one 
can find a Lagrangian radius at which the relaxation time 
(estimated in a manner to be described) is the same for all 
models at the same time. In a sense this can be interpreted as 
the radius which drives the expansion. In order to estimate 
the relaxation time at a given Lagrangian radius r we use 
eq.(2) with ru replaced by r. This does not take into account 
the fact that the mass within radius r varies with r, but since 
our purpose is to compare different models at the Lagrangian 
radius corresponding to the same fractional mass, this is not 
important. Fig. ^| shows the ratio of the relaxation times of 
the IK and 8K clusters at 4 different radii as a function of 
time. Throughout the post-collapse evolution, both clusters 
have similar relaxation times at the 10% radius, although 
there is some tendency towards larger radii towards the end. 
Similar results are obtained from a comparison of other N. 
The expansion of isolated clusters may therefore be thought 
of as being driven from around the 10% radius. The cluster 
cores will adjust themselves and their energy generation to 
deliver the necessary energy, while relaxation times at radii 
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Figure 3. Evolution of the lagrangian radii for clusters with dif- 
ferent initial particle numbers. Times are shifted such that all 
clusters go into core collapse at the same time. Low- TV clusters 
expand quicker after core collapse, and due to this faster expan- 
sion their relaxation times reach those of high-TV clusters. After 
10 5 TV-body time units, all clusters do therefore expand with the 
same rate. 



Figure 4. Fraction of bound mass as a function of TV-body time. 
Shown are clusters containing between TV = 128 to TV = 8192 
stars. Although low- TV clusters have initially smaller relaxation 
times and start losing mass earlier, curves for different TV meet 
each other at about 10 7 TV-body time units. Afterwards, high- 
TV models contain a smaller amount of mass at any given time, 
except for models with TV < 512. 



outward from the 10% radius are too large to influence the 
expansion much. 

Although it is clear from Fig. ^ that the post-collapse 
expansion cannot be strictly self-similar, the inner la- 
grangian radii up to the half-mass radius expand regular 
enough to be fitted by a relation like eq. 3. Fitting the evo- 
lution of different lagrangian radii and the mass-loss curve 
between 10 3 TV-body time units and 10 s TV-body time units 
(10 4 to 10 8 for TV ^ 4096) gives the ^-values shown in Ta- 
ble 2. For all TV, v drops with increasing distance from the 
cluster center, so that the core radii show the fastest ex- 
pansion. The reason for this will be discussed in the next 
section. 

Our values for the expansion of the half-mass radius 
agree fairly well with the values found by Goodman (1984) 
and Giersz & Heggie (1994b), despite the fact that different 
times were used by Giersz & Heggie (1994b) to obtain them 
and the scatter in v is fairly large for our high- TV runs. Un- 
fortunately no comparison values are available for the inner 
lagrangian radii. The best agreement with the (/-values de- 
rived from the mass-loss curve is achieved for the 10% and 
20% lagrangian radii. Since escapers are created mainly be- 
tween these radii (see section 3.4), a close correspondence in 
the behavior with time might be expected. 

Another interesting result of Table 2 is that the v val- 
ues derived from the mass-loss curve increase with TV, which 
means that high-TV clusters lose mass faster than clusters 
with smaller particle numbers. This is in contrast to a naive 
application of standard relaxation theory, in which it is as- 
sumed that v is independent of TV (eqs. 4). A large part 




R/R H .„ 



Figure 5. Density profiles for clusters starting with TV = 128 to 
TV = 8192 stars. Radii are given in units of the half-mass radius 
and densities are adjusted to be the same for all TV at the half- 
mass radius. The solid line shows a density profile proportional 
to p oc i? -3,29 . It gives a good fit to the halo profiles, especially 
for large TV. 
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Table 2. The parameter v determined from the mass-loss and 
the expansion rate of four different lagrangian radii as a function 
of the particle number N. Closest agreement with a self-similar 
expansion is achieved if the expansion of the 10% radius is com- 
pared with the mass-loss rate. 



N 


128 


256 


512 


1024 


2048 


4096 


8192 


vm 


0.10 


0.11 


0.11 


0.13 


0.13 


0.14 


0.14 


VrZ 


0.16 


0.16 


0.17 


0.16 


0.15 


0.17 


0.22 


1V10 


0.13 


0.11 


0.14 


0.14 


0.13 


0.11 


0.13 


V r 20 


0.08 


0.12 


0.13 


0.11 


0.08 


0.05 


0.09 


Vr50 


0.08 


0.09 


0.07 


0.02 


0.01 


0.06 


0.09 



of this discrepancy can be explained by the higher concen- 
tration of high-iV clusters, which breaks the homology and 
makes ki (eq. 3), and therefore v, increase with N. 

Plotting the bound mass as a function of time (Fig. 
shows that low-iV clusters start losing mass quicker due 
to their earlier core collapse, but that the increase in v is 
sufficient to reverse this after about 10 7 iV-body time units. 
Afterwards, high-iV runs contain a smaller fraction of their 
initial mass at a given JV-body time, except for runs with 
N < 512 which are affected by the increasing incompleteness 
towards the end. 

Finally, Fig. (and also Fig. ^|) shows that our 8k- 
model undergoes core oscillations up to about 10 4 iV-body 
time units, when the number of cluster stars drops below 
6700. The minimum required mass for core oscillations seems 
therefore to be of order 7.000, which agrees very well with 
results obtained by Goodman (1987) from his conducting 
gas sphere calculations. Only slightly larger values were ob- 
tained by Drukier et al. (1999) and Breeden et al. (1994), 
who found a lower limit of N ^ 8000. 

3.3 Density profile in post-collapse 

Fig. ^| depicts the relative sizes of the different radii for the 
2K and 8K models in greater detail. The time of the 8K 
model is rescaled by the initial relaxation time to match the 
relaxation time of a 2K model, and radii are scaled by the 
actual 20% radius. The build up of the halo is happening 
very slowly for both N and is complete only after about 10 5 
iV-body time units, corresponding to 10 3 initial relaxation 
times. By this time, the structure of the clusters outside 
the 20% radius is very similar for both N. Afterwards, the 
outermost radii in the 2K model start to decrease. This de- 
crease is caused by the recoil motion of the clusters due to 
escaping stars. Due to this recoil, weakly bound stars in the 
cluster halos gain positive energies relative to the cluster 
centers and escape. The recoil motion is larger and more 
erratic for low-TV clusters since a single star carries a larger 
fraction of the total cluster mass. Their halos will therefore 
be truncated at smaller radii. 

In contrast to the halo, the post-collapse profile in- 
side the half-mass radius is established already during core 
collapse and changes only little afterwards. The core radii 
of high-iV clusters are more concentrated relative to the 
20% radius than those of low-A r clusters. Giersz & Spurzem 
(1994) found that in low-iV models binaries have a higher 
chance to form and have to release less energy to prevent the 
core from contracting. They studied mainly the pre-collapse 
phase, however, given the similar density profile, the same 



i i i i i i i i i i 




0.1 r 



100 1000 10" 10 s 10 B 10 7 10" 10 s 10 10 10" 
T [NBODY] 



Figure 6. Ratio of the lagrangian radii relative to the 20% radius 
for clusters with N = 8192 (dashed line) and N = 2048 (solid line) 
stars. The time in the 8K model is rescaled by the relaxation time 
to that of the 2K model. Halos expand after core collapse until 
their build-up is reversed by the removal of weakly bound stars 
from the recoil motion of the center due to escaping stars. The 
cores are expanding slowly throughout the simulation because of 
increased binary activity in the center. 



conditions should prevail in the post-collapse phase. The 
larger core size of \ow-N clusters is therefore due to the 
stronger binary activity in their centers. Consequently, as 
the clusters evolve and lose mass, their core radii increase 
with time. 

Fig. |^ shows the density profiles at the time the clus- 
ter halos reach their maximum expansion relative to the 
20% radius. Radii are given in units of the half-mass radius 
and densities are adjusted such that they are the same at 
the half-mass radius. In the halo, the density profiles of all 
clusters follow power- law distributions p oc R~ a with slope 
a = 3.29. Deviations from this occur only in the outermost 
bins, where small- N clusters have steeper profiles because 
of the recoil motion from escaping stars, which effectively 
truncate the clusters at a certain radius. 

Our value for a is not far from the theoretical prediction 
by Spitzer & Shapiro (1972), who found a = 3.5 for the 
halo density of an isolated star cluster. Similar slopes were 
also found by Lightman & Shapiro (1978) and by Takahashi 
(1996) in his Fokker-Planck simulations of anisotropic post- 
collapse clusters. 

The transition between the halo and the inner parts 
happens at around R = 0.4 half-mass radii. Inward from 
there, clusters have shallow density profiles which become 
steeper for increasing N. Inagaki & Lynden-Bell (1983) and 
Takahashi (1995) found that the inner profiles can be fitted 
by power-law distributions with slopes of a — 2.0 and a = 
2.23 respectively. Our profiles do not strictly follow power- 
laws, but can be fitted by them over a restricted range in R. 
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Figure 7. Cumulative mass distribution at half-mass time for 
clusters starting with different density profiles. Clusters starting 
with Plummer profiles are shown by a solid line, clusters starting 
with King Wo = 3.0 profiles by a dashed line. The distributions 
at half-mass time agree approximately with each other since the 
cluster expansion has erased the memory of the starting condition. 

For the 4K and 8K models, an increase of p ~ R~ 233 gives 
a good fit between 0.02 and 0.3 half-mass radii, which is not 
too different from the result of Takahashi (1995). Given the 
relatively small number of stars in our simulations, and the 
fact that the profiles have not reached a limiting profile even 
for our largest N, it might be possible that a single power- 
law profile is reached in the centers of large- N clusters. 

3.4 Evolution along a single sequence of models ? 

In the expansion phase after core collapse, the space that 
the clusters initially occupied becomes a vanishing fraction 
of their actual volume. One might therefore ask if the density 
profile after a large enough time has passed still depends on 
the initial profile. In order to test this assumption, we made 
a set of 8 comparison runs which started from King Wo = 3.0 
profiles and had N = 1024 stars initially. 

Fig. compares the cumulative mass distribution of the 
King-models with the TV = 1024 Plummer-model, after both 
have lost half their initial mass (note that the clusters have 
expanded by about a factor of 100 by this time). A cumu- 
lative distribution has the advantage that it better shows 
the behavior of the bulk of the stars. It can be seen that the 
mass profiles of both clusters agree very well with each other. 
Similar comparisons at other times give the same agreement, 
provided the time is much larger than the core collapse time. 
The density distribution of isolated clusters long after core 
collapse becomes therefore independent of their initial den- 
sity profile. 

Similarly one might ask whether the density profile or 
any other cluster parameter at a given particle number still 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

M(<R)/M Tol M(<R)/M Tol 



Figure 8. Distribution of density and anisotropy {3 of clusters 
starting with different initial particle numbers by the time 512 
stars are bound. Radii are scaled by the half-mass radius. /3 is cal- 
culated as explained in section 3.7. Both the mass and anisotropy 
profiles arc independent of the initial particle number TVq . 
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Figure 9. Evolution of the number of bound stars for different 
clusters. Radii are rescaled such that the size of high- TV clusters 
when they have 45% bound stars is equal to the size of low-TV 

clusters with 90% bound stars. TV-body time units are scaled by 

3/2 

r h . After this rescaling, the bound mass evolves similar in both 
cases, showing that the evolution is independent of the initial 
particle number. 

depends on the initial number of cluster stars. Fig. ^ depicts 
the density and anisotropy profile as a function of enclosed 
mass by the time 512 stars remain bound, for clusters start- 
ing between TV = 1024 to TV = 4096 stars. The density 
and anisotropy profiles agree well with each other for any 
radius, and similar plots at other times also give excellent 
agreement between the runs. 

We conclude that the cluster structure becomes inde- 
pendent of the initial particle number and density distribu- 
tion after enough time has passed and is solely a function 
of the actual particle number. Clusters in late post-collapse 
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can therefore be characterised by only two parameters, as 
e.g. the current particle number and the half-mass radius. 
This would imply that isolated clusters evolve along a single 
sequence of profiles and that the evolution of clusters which 
started from different initial configurations is similar. 

To check this assumption, we have a look at the evo- 
lution of bound mass with time. If clusters evolve along a 
single sequence, the evolution of low- TV clusters should be 
similar to that of high- TV clusters when they have reached 
the same number of stars. Fig. ^compares the mass-loss rate 
of small- TV clusters with clusters that had initially twice as 
many stars. High-iV runs are rescaled such that the clus- 
ter sizes when 45% of the stars remain bound are equal to 
the sizes of the low- TV" clusters when 90% of their stars are 
bound. The time is also shifted such that clusters reach 90% 
and 45% bound stars at the same time. It can be seen that 
we obtain a very good fit for the evolution of bound mass 
with time. The fit is best when many simulations were made 
as e.g. in cases with TV = 512 and TV = 1024. Deviations oc- 
cur, if at all, only at the end of the runs when the particle 
numbers are small. We conclude that isolated clusters evolve 
along a single sequence of models. 
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3.5 Escape of stars 

It was shown by Henon (1960) that stars can escape from 
isolated clusters only due to close encounters with other 
stars or binaries. If close encounters are the dominant es- 
cape mechanism, most stars should escape from the inner 
cluster regions, where the density of stars is highest and 
close encounters between stars should happen most often, 
though the potential well is deepest. 

In our simulations, we checked the energy of each star 
when its regular force in NBODY6 was calculated. In order 
to do this, potential energies were calculated with respect 
to all stars and the kinetic energy was calculated relative to 
the cluster velocity from the time of the last data output. 
The time and radius where stars acquired positive energies 
were noted and we will refer to this radius as the 'scattering 
distance' Rscatt- If stars still had positive energies after their 
distance from the cluster center had increased by a factor of 
two compared to the scattering distance, we assumed that 
the star was successfully scattered out of the cluster and 
measured its excess energy above the E = energy thresh- 
old. Stars that fell below E — before they reached twice 
the scattering distance were not counted as escapers. 

Fig. [To] shows, as an example for a low- TV cluster, the 
distribution of energies against scattering radius for the 
TV = 512 star models. The energies are expressed in terms 
of the one-dimensional velocity dispersion of all bound stars 
lkT = 2/3 <Estar> and the distances are measured in 
units of the 10% radius. One can notice three different con- 
centrations in this diagram. Stars escaping due to single-star 
encounters have typical energies of a few tenths of lkT and 
acquire the escape energies at distances around the 10% la- 
grangian radius. Single-star escapers coming from distances 
closer to the cluster center have higher energies since stars 
at smaller radii move with larger velocities and can therefore 
also create higher energy escapers. The most distant single 
star escapers come from around R = 10 R\o% , which corre- 
spond roughly to the half-mass radius. Stars escaping due 
to single-star encounters form the majority of escapers. 



Figure 10. Distribution of escape energies against radius where 
stars become unbound for the N = 512 star clusters. Plotted 
are all stars that escape before the half-mass time. The escape 
energies are measured after the stars have reached twice the scat- 
tering distance. Most stars escape due to encounters with single 
stars around the 10% lagrangian radius. Escapers with high en- 
ergies come from encounters with binaries in the cluster cores. 
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Figure 11. Same as Fig. |lG| for the N = 8192 cluster. The distri- 
bution is similar to the N = 512 case, except that escapers from 
binary encounters are coming from much smaller distances. 
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Figure 12. Distribution of escape energies as a function of scaled time for clusters with N = 512 stars. Time is scaled according to eq. 
1. The dotted line shows the distribution of all stars, the solid line shows the distribution after stars that escape due to the movement of 
the cluster have been subtracted. Escapers due to single star encounters follow roughly a gaussian distribution. The dashed line shows 
the remaining distribution after this gaussian has been subtracted from the high-energy escapers. Binary induced escape sets in only 
after core collapse and also follows roughly a gaussian distribution in energy. 



Stars which escape because of binary-single star inter- 
actions are created only in the cluster cores and have typical 
energies between 2 and 100 kT. Since the energy of the es- 
caper is coming from the internal energy of the binary, there 
is no dependence of the escape energy on the scattering ra- 
dius. Although smaller in number, escapers due to binary 



encounters carry away more than 95% of the total energy of 
escapers. 

A third group of stars is created in the halos at distances 
larger than 10 half-mass radii. These stars escape since the 
cluster center moves erratically as a result of the recoil from 
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high energy encounters in the core, which will unbind loosely 
bound stars. 

A similar plot for the 8K model (Fig. [H]) shows that 
the energy of the single star escapers changes very little with 
particle number. In addition, the mean scattering distance 
remains the same. Escapers due to binaries are now well 
separated from the single star escapers. Their mean distance 
from the center has decreased by almost a factor of 10 due 
to the more concentrated core of the high- TV model. 

Fig. [li] shows the distribution of escape energies as a 
function of scaled time (calculated according to eq. 1) for 
TV = 512. Escape of high-energy stars sets in only after core 
collapse (happening at around T = 150 TV-body time units). 
In agreement with Giersz & Heggie (1994a), we find that be- 
fore core collapse, stars escape due to single star encounters 
only. The typical escape energies of single-star escapers show 
no variation with time and lie between 0.04 and lkT. After 
stars escaping due to the movement of the cluster are taken 
out, the energy distribution can be fitted by a sum of two 
gaussians, which correspond to the two escape mechanisms. 
The majority of stars, typically about 80% to 85%, escape 
due to single-star encounters. 

Fig. his shows the escape energies and escape times for 
stars in the 8K model. The energy distribution of stars es- 
caping due to single-star and binary encounters is similar to 
the TV = 512 case. Escapers by binary encounters tend to 
have slightly higher energies than in the TV = 512 star case. 
The fraction of binary escapers is slightly smaller than in the 
TV = 512 case because the fraction of binaries compared to 
the number of cluster stars has decreased. The most striking 
difference is that in the 8K model, stars from binary encoun- 
ters do not have a smooth distribution in time, but escape 
in several distinct phases. 

Fig. [l4| shows that these phases are linked to the core 
oscillations of the 8K model. Stars escape by binary encoun- 
ters only when the cluster is in core collapse since the core 
density is much lower in an expansion phase where close 
encounters between binaries and single stars are unlikely. 
Stars escaping because of single-star encounters show much 
less variation with time as they are typically created at 
larger radii that change less due to the oscillations. Fig. [l^ 
also shows that the rate at which single-star escapers are 
created is increasing during and after core collapse which 
could be the result of the increasing anisotropy. As a con- 
sequence, more loosely bound halo stars will travel through 
high-density regions where they have a high chance of being 
scattered to escape energies. In addition, the overall change 
in the density distribution can also increase the escape rate. 

3.6 Binaries 

Fig. [IB] shows the number of bound binaries as a function 
of time. We assume that the number of binaries is equal to 
the number of regularised pairs in NBODY6, which should 
be the case except for transient encounters between cluster 
stars which are also regularised if they are close enough. 

No binaries are present at the start of our simulations 
and they are created only during and after core collapse. 
The maximum number of bound binaries is reached after 
about 10 core collapse times. Throughout the post-collapse 
expansion, a few binaries are sufficient to drive the expansion 
of the clusters, even for clusters with several thousand stars. 
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Figure 14. Evolution of the escape energies (bottom panel), core 
radius (middle) and binary number (top) in the 8K model. Al- 
though binaries are present all the time, they are efficient in cre- 
ating escapers only when the cluster is in core collapse. 

The number of binaries increases with the number of 
cluster stars, while their fraction decreases. If the number of 
binaries is taken at the maximum value, the increase can be 
fitted approximately as Nsin oc TV 0,3 . This is comparable to 
the increase given by Giersz & Heggie (1994b), who found 
Nsin oc TV 0,2 after a few core collapse times had passed. 
Given the considerable scatter in our data and the fact that 
the binary number in our models is still rising after a few 
core collapse times, the agreement is satisfactory. 

Comparison of the members of binary systems at suc- 
cessive times shows that the vast majority of them are not 
transient encounters but bound systems, especially for clus- 
ters with high particle numbers. In addition, the number of 
binaries in the core is increasing with the particle number 
in a similar way than the total binary number. Hence the 
number of active binaries needed to support the expansion of 
the clusters seems to rise with TV, in contrast to Goodman's 
(1984) results from gaseous model calculations. 

For later times, the binary fraction starts to drop again 
as the clusters contain fewer and fewer stars. This effect can 
best be seen for the high- TV runs. Runs with small particle 
numbers show almost no decrease, probably because at least 
one binary has to present in the cluster core to drive the 
evolution. If this star is ejected, a new binary will be formed 
quickly, so a value of one forms a lower limit for the number 
of binaries (cf. Fig. 15). 

The radial distribution of bound binaries is shown in 
Fig. |l| Although more concentrated then single-stars, a sig- 
nificant fraction of binaries can be found outside the core. 
Only one quarter of all binaries is located within the 3% la- 
grangian radius. Binaries outside the core are mostly ejected 
by three and four body interactions and move on very elon- 
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Figure 13. Same as Fig. [12| for N = 8192 stars. The distribution of escapers is very similar to the N = 512 case. However, stars that 
escape due to encounters with binaries do not escape continuously, but in several distinct phases which are related to the core oscillations 
of the cluster. 



gated orbits. Subsequent passages through the cluster center 
will remove them from the cluster. Given the smaller density 
of stars and binaries in the halo, binary activity should be 
confined to the cluster centers. 



3.7 Anisotropy 

Fig. jnj shows the evolution of the mean anisotropy for clus- 
ters with N — 4096 stars. To calculate the anisotropy, only 
bound stars are used. Following Binney & Tremaine (1987), 
we define the anisotropy parameter as f3 = 1 — (<rf/2cr^), 
which is half as large as the A-parameter used by Giersz 
& Spurzem (1994) and Takahashi (1996). The anisotropy 
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Figure 15. Number of bound, regularised binaries as a function 
of time. The binary number reaches a maximum after about 10 
core collapse times. It is larger for high-TV clusters and decreases 
as the clusters lose stars. 
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Figure 16. Radial distribution of binaries and single-stars for 
clusters containing TV = 512 and TV = 2048 stars respectively. 
Binaries (dashed) are more concentrated than single stars (solid 
line). In addition, binaries in high- TV clusters can be found at 
smaller radii. 



rises from the start of the simulations and reaches a maxi- 
mum after about 10 6 iV-body time units, corresponding to 
1.5- 10 4 initial relaxation times. We find that the same num- 
ber of initial relaxation times is required for other TV. We 
also find a close agreement between the time when the clus- 
ters reach maximum anisotropy and the time when their 
halos reach maximum expansion relative to the inner radii. 
At this time, the outermost shells have become nearly 100% 
radially anisotropic. Clusters remain isotropic inside their 
20% radius throughout the simulation, which agrees very 
well with results found by Takahashi (1996) from Fokker- 
Planck simulations. 

One can see that the anisotropy is decreasing for very 
late times, which is the result of the increasing removal of 
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Figure 17. Evolution of the anisotropy with time for 4 differ- 
ent radii for clusters with TV = 4096 stars. Clusters are strongly 
anisotropic in their outer parts but remain isotropic near their 
centers. The mean anisotropy increases up to 10 6 TV-body time 
units and decreases slowly afterwards. A similar behavior is found 
for other particle numbers. 

weakly bound stars as the particle numbers drop. Stars with 
energies near E = are predominantly on near radial orbits 
when closer to the cluster center. The removal of such stars 
must therefore decrease the anisotropy at intermediate radii. 

We can compare the anisotropy profile for clusters with 
different particle numbers by binning data from several 
snapshots around the time when the clusters reach maxi- 
mum anisotropy. Fig. |l8| shows the resulting profiles as a 
function of enclosed mass. The anisotropy rises monoton- 
ically from the core, which is isotropic for all TV, towards 
the halo. For large TV, the differences become smaller and it 
might be possible that the profiles are converging towards a 
maximum profile at the high- TV end. Unfortunately it is dif- 
ficult to say if this is really the case, due to the small number 
of simulations made (only one for TV = 8192) and the result- 
ing uncertainties in the profiles. If so, it would mean that 
there is an universal density and velocity profile, which only 
clusters with large enough TV can reach, since the evolution 
of low- TV models towards this profile is interrupted by the 
removal of their outer halo. 

It was already shown that the cluster expansion is 
driven by encounters of stars near the 10% lagrangian ra- 
dius (Fig. ^|). Hence, the 10% lagrangian radius should also 
play an important role for the anisotropy profile. Fig. |l^ 
depicts the anisotropy profile for different particle numbers, 
calculated by dividing the radii of each model by its 10% 
lagrangian radius. It can be seen that the profiles are nearly 
the same for all TV, showing again that the 10% radius is the 
radius where halo stars are created. We found that no other 
radius leads to such an agreement in the profiles. 
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Figure 18. Maximum anisotropy in dependence of the initial 
number of cluster stars. The anisotropy is plotted as a function of 
enclosed mass. Shown are curves from N = 256 (bottom) to N = 
8192 (top). The anisotropy in the halo increases with increasing 
particle number. 
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Figure 19. Maximum anisotropy as a function of distance to 
the cluster center. If distances are scaled by the 10% lagrangian 
radius, the anisotropy profile is independent of the particle num- 
ber. Most stars should therefore be scattered into the halo from 
around this radius. 



Henon (1973) found that spherical stellar systems with 
radially anisotropic distribution functions can be vulnerable 
to radial orbit instabilities. His results were later confirmed 
and extended by Barnes, Hut & Goodman (1986) and De- 
jonghe & Merritt (1988). Dejonghe & Merritt (1988) found 
that models based on Plummer's density law become un- 
stable if the global anisotropy exceeds a critical value of 
2T r /Tt « 2.0. The anisotropies in our highest- N models are 
close to this value, and so it might be interesting to look for 
signs of such instabilities in our runs. 

We calculated for all models with N ^ 1024 and all 
times when data was stored the moment of inertia tensor of 
the stellar distribution and its eigenvectors and eigenvalues. 
We studied only radial shells between the 20% and 80% la- 
grangian radius, since the velocity distribution is isotropic 
for smaller radii. Since previous work indicates that radially 
anisotropic stellar systems develop bar-like mass distribu- 
tions, we first checked for stellar bars in our models. 

Fig. [20] shows for the 8K run the angular distribution 
of the eigenvectors which belong to the smallest moment of 
inertia. It can be seen that their distribution is more or less 
random. This is not what one would expect in the presence of 
radial orbit instabilities, since bar-like configurations should 
have fixed axes in space, at least for a certain amount of 
time. In a second test we compared the ratio of the smallest 
eigenvalue to the sum of the other two. If a cluster collapses 
into a bar-like configuration, one would expect to see a sud- 
den drop of this ratio. However, no such drop could be seen 
anywhere in the run. Hence no signs of instabilities could be 
found in the 8K run and similar results were obtained for 
runs with smaller particle numbers. We also found no signs 
for disc-like instabilities in any of our runs. 

There are several possible explanations for the absence 
of radial orbit instabilities: Our distributions might not be 
anisotropic enough to become unstable, or relaxation pro- 
cesses suppress them since they tend to make distribu- 
tions spherical. Most important for halo stars are passages 
through the inner cluster parts, where encounters with other 
stars can significantly change their orbits. As was shown 
earlier (Fig. 3), the relaxation times in the center change 
only slowly with the particle number since high- TV clusters 
are more concentrated. The onset of radial orbit instabilities 
might therefore be suppressed by relaxation processes in the 
cluster center, no matter how large TV becomes. 

A third possibility is that the halo is constantly renewed 
and extended by stars ejected from the cluster center. Since 
such stars will be scattered uniformly in all directions, they 
can also prevent the formation of structures in the halo. 



4 CONCLUSIONS 

We have performed a set of TV-body simulations of isolated 
single-mass clusters starting from Plummer profiles. Our 
main focus was on the post-collapse evolution, and for the 
first time we could follow this evolution until nearly com- 
plete disintegration of the clusters, thereby extending the 
parameter space of collisional stellar dynamics significantly. 

We found that after a sufficient number of relaxation 
times has passed, the structure of the clusters becomes in- 
dependent of the initial density profile and particle number, 
and can be characterized by just two parameters, as e.g. the 
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Figure 20. Angular distribution of the eigenvectors that belong 
to the smallest eigenvalue for the 8K run. Circles show the pro- 
jection into the x-y plane, open triangles into the x-z plane and 
filled triangles into the y-z plane. The solid line shows the ex- 
pected value for a uniform distribution. The distribution of the 
data from the TV-body run is compatible with a random distribu- 
tion. 

actual particle number and the half-mass radius. For very 
large N, the cluster structure might become independent 
of the actual particle number too, although our simulations 
could not reach large enough N to test this assumption. We 
found that isolated clusters evolve along a single sequence 
of models and runs starting from different particle numbers 
can be stitched together. This can be relevant for the post- 
collapse evolution of galactic globular clusters which have 
half-mass radii much smaller than their tidal radii, since 
these systems would be nearly isolated and might have pro- 
files similar to our models. 

We found that low- TV clusters have larger radii after core 
collapse, so that all clusters have similar relaxation times 
during the post-collapse expansion phase. As a consequence, 
the mass-loss rates depend only weakly on the particle num- 
ber. Due to their higher concentration, high-A^ clusters loose 
mass even faster than low-A^ models. 

Stars in isolated clusters are ejected almost entirely 
from within the half-mass radius, with encounters between 
single stars being the dominant escape mechanism. We could 
establish a clear correlation between the excess energy of es- 
capers and the position where the escapers are created. Such 
a relation could also exist for clusters in tidal fields, although 
the long time necessary for escape complicates the escape of 
stars in this case (Fukushige & Heggie 2000, Baumgardt 
2001). 

For all clusters, only a few binaries are present and drive 
the expansion. Their fraction decreases with increasing par- 
ticle number, while the total number increases roughly as 
Nsin ~ A' ' 3 . About 15% of the stars are ejected due to 



encounters with binaries and these carry away a large frac- 
tion of the cluster energy. Binaries are strongly concentrated 
towards the cluster cores, and the ejection by binaries also 
happens only near the cluster centers. For clusters show- 
ing core oscillations, binary induced escape is efficient only 
in the contraction phases of the core, while the escape due 
to single-star encounters fluctuates much less with the core 
oscillations. 

A radially anisotropic velocity distribution is created 
in the cluster halos, mainly as a result of two-body en- 
counters near the 10% lagrangian radius. High-A clusters 
have smaller core sizes and show therefore larger overall 
anisotropies. It seems conceivable that isolated clusters are 
vulnerable to radial orbit instabilities for sufficiently large 
N, but no indication for such instabilities could be found. 
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